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Abstract 

In this paper we will address the problem of developing mathematical models for 
the numerical simulation of the human circulatory system. In particular, we will focus 
our attention on the problem of haemodynamics in large human arteries. 
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1. Introduction 

The simulation of not only the physiological functioning of the blood circulatory 
system, but also of specific pathological circumstances is of utmost importance since car- 
diovascular diseases represent the leading cause of death in developed countries, with a 
tremendous medical, social and economic impact. 

In the cardiovascular system, altered flow conditions, such as flow separation, flow 
reversal, low and oscillatory shear stress areas, are recognised as important factors in the 
development of arterial diseases. A detailed understanding of the local haemodynamics, 
the effect of vascular wall modification on flow patterns and its long-term adaptation to 
surgical procedures can have useful clinical applications. Some of these phenomena are 
not well understood, making it difficult to foresee short and long term evolution of the 
disease and the planning of the therapeutic approach. In this context, the mathematical 
models and numerical simulations can play a crucial role. 

Blood flow interacts both mechanically and chemically with the vessel walls. The me- 
chanical coupling requires algorithms that correctly describe the energy transfer between 
the fluid (typically modelled by the Navier-Stokes equations) and the structure of the vessel 
wall. On the other hand, the flow equations can be coupled with appropriate models that 
describe the wall absorption of bio-chemicals (e.g. oxygen, lipids, drugs, etc.) and of their 
transport, diffusion and kinetics. Numerical simulations of this type may help to understand 
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the modifications in bio-chemical exchanges due to an alteration of the flow field caused, 
for instance, by a stenosis (i.e. a localised narrowing of a vessel lumen, normally due to fat 
accumulation). 

The simulation of large and medium-size arteries is now sufficiently advanced so to 
envisage the applications of computer models to medical research and, in a medium range, 
to everyday medical practise. For instance, simulating the flow in a coronary by-pass may 
help understanding the extent at which its geometry influences the flow and in turn the 
post-surgery evolution. Also the study of the effects of a vascular prosthesis as well as the 
study of artificial valve implants are areas which could benefit from a sufficiently accurate 
simulation of blood flow field. 

In this paper we review the principal mathematical steps behind the derivation of the 
coupled fluid-structure equations which model the blood flow motion in large and medium- 
sized arteries. Then we mention the way geometrical multiscale models, that combine 
mathematical models set up in different spatial dimensions, can be conveniently used to 
simulate the whole circulatory system. 



In this section we will treat the situation arising when the flow in a vessel interacts 
mechanically with the wall structure. This aspect is particularly relevant for blood flow 
in large arteries, where the vessel wall radius may vary up to 10% because of the forces 
exerted by the flowing blood stream. 

We will first illustrate a framework for the Navier-Stokes equations in a moving do- 
main which is particularly convenient for the analysis and for the set up of numerical solu- 
tion methods. 

2.1. The Arbitrary Lagrangian Eulerian (ALE) formulation of the 
Navier-Stokes equations 

Navier-Stokes equations are usually derived according to the Eulerian approach where 
the independent spatial variables are the coordinates of a fixed Eulerian system. When 
considering the flow inside a portion of a compliant artery, we have to compute the flow 
solution in a computational domain il t varying with time. 



Figure 1 : A simple model of a section of an artery. The vessel wall Tf 
is moving. The location along the z axis of and are fixed. 

The boundary of Q t may in general be subdivided into two parts. The first part coin- 
cides with the physical fluid boundary, i.e. the vessel wall T™ (see Fig.Q, which is moving 



2. The coupled fluid-structure problem 



r, 
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under the effect of the flow field. The other part of dQ, t corresponds to artificial boundaries 
which delimit the region of interest from the remaining part of the cardiovascular system. 

The "artificial" boundaries are the inlet and outlet (or, using the medical terminology, 
the proximal and distal) sections, here indicated by T" and Tf*, respectively. The location 
of these boundaries is fixed a priori. More precisely, T" and Tf^ may change with time 
because of the displacement of Ff, however they remain planar and their position along 
the vessel axis is kept fixed. In this case the Eulerian approach becomes impractical. 

A possible alternative would be to use the Lagrangian approach, where we identify 
the computational domain on a reference configuration Q , e.g. that at the initial time 
t = 0, and the corresponding domain in the current configuration will be provided by the 
Lagrangian mapping 

Q t = Qc t =C t (Q ), t>0, (2.1) 

which describes the motion of a material particle and whose time derivative is the fluid 
velocity. Since the fluid velocity at the wall is equal to the wall velocity, the Lagrangian 
mapping effectively maps to the correct wall position at each time t. However, the 
artificial boundaries in the reference configuration, say and T^, will now be transported 
along the fluid trajectories. This is unacceptable, particularly for a relatively large time 
interval as Cl t rapidly becomes highly distorted. 

A more convenient situation is the one when, even if the wall is moving, one keeps 
the inlet and outlet boundaries at the same spatial location along the vessel axis. With 
that purpose, we introduce the Arbitrary Lagrangian Eulerian (ALE) mapping At : Oo — ► 
il^4 t , Y — > x(t, Y) = At(Y), which provides the spatial coordinates (i, x) in terms 
of the so-called ALE coordinates (t, Y), with the basic requirement that At retrieves, at 
each time t > 0, the desired computational domain, £l t = £1.4 t = ^4 t (fi ), t > 0. 

The ALE mapping should be continuous and bijective in il - Once given, we may 
define the domain velocity field as 

w=^oA-\ (2.2) 
where the composition operator applies only to the spatial coordinates. The ALE time 

D± 
Dt 



derivative of a function / : 7 x fi ( -> R, which we denote by -7^7-/, is defined as 



/":/xfi t ->R, f = -±oA7 1 . (2.3) 

This definition is readily extended to vector valued functions. The ALE derivative is related 
to the Eulerian (partial) time derivative by the relation 

D A df 

w / = i +w - v/ ' (2 - 4) 

where the gradient is made with respect to the x-coordinates. 

The Navier-Stokes equations may be formulated in order to put into evidence the 
ALE time derivative, obtaining 

+ [(u - w) • V] u + Vp - divT(u) = f, 

in n t , t>0, (2.5) 

divu = 0, 
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where T is the Cauchy stress tensor, which for Newtonian fluid is given by T = i/(Vu + 
Vu T ), being v the blood kinematic viscosity. 

When considering small vessels, accounting for non-Newtonian behaviour of blood 
becomes crucial. In that case the functional dependence of T on u becomes more complex, 
see for instance 1 18 1. 



P™sA ~ divcr (v, si) =f, inftf, t>0, 



2.2. The structure model 

The vascular wall has a very complex nature and devising an accurate model for its 
mechanical behaviour is rather difficult. Its structure is indeed formed by many layers with 
different mechanical characteristics Ri ll II . which are usually in a pre-stressed state. More- 
over, experimental results obtained by specimens are only partially significant. Indeed, the 
vascular wall is a living tissue with the presence of muscular cells which contribute to its 
mechanical behaviour. It may then be expected that the dead tissue used in the laboratory 
will have different mechanical characteristics than the living one. Moreover, the arterial 
mechanics depend also on the type of the surrounding tissues, an aspect almost impossible 
to reproduce in a laboratory. As we have already pointed out, the displacements cannot be 
considered small (at least in large arteries where the radius may vary up to a few percent 
during the systolic phase). Consequently, an appropriate model for the structure displace- 
ment tj reads 

— -divTfo,- 

where fif indicates the current configuration and <x is the Cauchy stress tensor. The latter 
may depend on the structure velocity because of viscoelasticity. A full Lagrange formula- 
tion for the structure on a fixed reference configuration fig may be obtained by the usual 
Lagrange and Piola transformation (see, e.g., 1 2 1). A general framework to derive constitu- 
tive equations for arterial walls is reported in II II . 

It is the role of mathematical modelling to find reasonable simplifying assumptions 
by which major physical characteristics remain present, yet the problem becomes compu- 
tationally actractive. In particular, a simpler model may be obtained by considering only 
displacements in the radial direction and a cylindrical geometry for the vessel. Furthermore 
if we neglect the geometrical non-linearities (which correspond to assume small displace- 
ments), as well as the variations along the radial directions (small thickness assumption) 
we obtain the following "generalised string model" 1 16 14 1 for the evolution of the radial 
displacement r) = R — Rq, 

where a > and b > are parameters linked to the vessel geometry and mechanical 
characteristics, c > is a viscoelastic parameter and H is a forcing term which depends on 
the action of the fluid, as we will see in ( 12.71 . More details as well as the derivation of this 
model are found in the cited references. 

Here Tq is the reference configuration for the structure 

Tg = {(r, 9,z):r = R (z), 9 e [0, 2tt), z e [0, L]}, 
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where L indicates the length of the arterial element under consideration. In our cylindrical 
coordinate system (r, 0, z), the z coordinate is aligned along the vessel axes and a plane 
z = z (= constant) defines an axial section. 



2.3. Coupling with the structure model 

We now study the properties of the coupled fluid-structure problem, using for the 
structure the generalised string model ( 12. 6> . We will take n always to be the outwardly 
vector normal to the fluid domain boundary. Furthermore we define g as the metric function 
so that the elemental surface measure da on Tf is related to the corresponding measure dao 
on Tq by da — g dao. 

We will then address the following problem: For all t > 0, find u, p , rj such that 

u, p satisfy problem ( 12. 51 . 
rj satisfies problem ( 12. 6> . 



uo At = ^e r , on r™, (2 ' 7) 
ff = 5-^-[(p-po)n-T(u)-n]-e r on 

Here, po is the pressure acting at the exterior of the vessel, e r is the radial unit vector, p w 
and p are the wall and fluid densities, respectively, while ho is the wall thickness. The 
system is complemented by appropriate boundary and initial conditions. 

We may then recognise the sources of the coupling between the fluid and the structure 
models, which are twofold. In view of a possible iterative solution strategy, the fluid solu- 
tion provides the value of H, which is function of the fluid stresses at the wall. On the other 
hand, the movement of the vessel wall modifies the geometry on which the fluid equations 
must be solved, besides providing Dirichlet boundary conditions for the fluid velocity in 
correspondence to the vessel wall. 

Remark 2.1. We may note that the non-linear convective term in the Navier-Stokes equa- 
tions is crucial to obtain the well-posedness of the coupled problem, because it generates 
a boundary term which compensates that coming from the treatment of the acceleration 
term. These two contributions are indeed only present in the case of a moving boundary. 
Seed and Q. 



2.4. Numerical solution of the coupled fluid-structure problem 

In this section we describe an algorithm that at each time-level allows the decoupling 
of the sub-problem related to the fluid from that related to the vessel wall. As usual, t k , 
k = 0, 1, . . . denotes the k-th discrete time level; At > is the time-step, while v k is the 
approximation of the function (either scalar or vector) v at time t k . 

The numerical solution of the fluid-structure interaction problem ( 12. 7> will be car- 
ried out by constructing a proper finite element approximation of each sub-problem. In 
particular, for the fluid we need to devise a finite element formulation suitable for mov- 
ing domains (or, more precisely, moving grids). In this respect, the ALE formulation will 
provide an appropriate framework. 
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To better illustrate the situation we refer to Fig. |2] where we have drawn a 2D fluid 
structure interaction problem (only the upper portion of the vessel is reported). For the sake 
of simplicity we have considered only a two-dimensional fluid structure problem, yet the 
algorithm here presented may be readily extended to more complex situations and three 
dimensional problems. 




L z 



Figure 2: A simple fluid-structure interaction problem. On the left the 
domain definition and on the right the discretized vessel wall corre- 
sponding to a possible value of r/h- 




Figure 3: The triangulation used for the fluid problem at each time t is 
the the image through a map At of a mesh constructed on Slo- 

The structure on Tq will be discretized by means of a grid T t s h and employing piece- 
wise linear continuous (PI) finite elements to represent the approximate vessel wall dis- 
placement rjh- The position at time t of the discretized vessel wall boundary, correspond- 
ing to the discrete displacement field is indicated by Tf h . Consequently, the fluid 
domain will be represented at every time by a polygon, which we indicate by Clt,h- Its 
triangulation T/ h will be constructed as the image by an appropriate ALE mapping At of 
a triangulation T^ h of flo> as shown in Fig. [3] Correspondingly, il t .h = At&o,h< where 

fio,7i is the discretisation of fio induced by TJ h . The trace of Tj h on Tq will coincide with 
the T t s h of the vessel wall, i.e. we here consider geometrically conforming finite elements 
between the fluid and the structure. The possibility of using geometrically non-conforming 
finite elements has been investigated in 1101 . 

Having at disposal the discrete displacement field at t — t k+1 and thus the 
position of the domain boundary dVL t h+i h , the set up of a map Afh+i such that A t k+i (T , h ) 
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is an acceptable finite element mesh for the fluid domain is not a simple task. However, if 
we can assume that il t ,h is convex for all t and that the displacements are relatively small, 
the following technique, known as harmonic extension, may well serve the purpose. If X/j 
indicates the PI finite element vector space associated to T/ h and g/j : 0Qo,h — > dil t k+i h 
is the function describing the fluid domain boundary, we build the map by seeking G 
such that 

/ Vy h : Vz h = Vz,, G X°, y h = Sh , on dQ ,h, (2.8) 
Jn 

and then setting A t k+i(Y) — y/i(Y), VY G f2o,ft. A more general discussion on the 
construction of the ALE mapping may be found in (5)11 21 as well as in |9|. 
Remark 2.2. Adopting PI elements for the construction of the ALE map ensures that the 
triangles of T^ are mapped into triangles, thus T^ t is a valid triangulation, under the 
requirement of invertibility of the map (which is assured if the domain is convex and the 
wall displacements are small). 

As for the time evolution, we may adopt a linear time variation within each time slab 
[t k ,t k+1 ] by setting 

f _ fk f _ j.fe+1 

A = —^j-Atk+x - — ^ — Atk, te[t k ,t k ^}. 

Then, the corresponding domain velocity w/, will be constant on each time slab. 



2.4.1. The iterative algorithm 

We are now in the position of describing an iteration algorithm for the solution of the 
coupled problem. As usual, we assume that all quantities are available at t = t k , k > 0, 
provided either by previous calculations or by the initial data and we wish to advance 
to the new time step t k+1 . For ease of notation we here omit the subscript h, with the 
understanding that we are referring exclusively to finite element quantities. 

The algorithm requires to choose a tolerance r > 0, which is used to test the con- 
vergence of the procedure, and a relaxation parameter < 8 < 1, In what follows, the 
subscript j > denotes the sub-iteration counter. 

The algorithm reads: 

Al Extrapolate the vessel wall structure displacements and velocity: 

V(o) =r ) + Atr l i »7( ) =V ■ 

A2 Set j = 0. 

A2. 1 By using rff) 1 compute the new grid for the fluid domain Q t and the ALE map 

by solving the harmonic extension problem (12.81 . 
A2.2 Approximate the Navier-Stokes problem to compute and p^ 1 -^, using 

as velocity on the wall boundary the one calculated from 

A2.3 Approximate the structure problem to compute "q k+1 and rj k+1 using 

and P(j^x) to recover tne forcing term H. 
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A2.4 Unless ||r/* 



,fe+i 



= erfa 1 + (1 - ^ = ^g 1 + (i - 0)>7* fc 



fe+i 



and j <— j + 1. Then return to step|2al 



A3 Set 



fe+i 



= V* 



fe+i 




u 



fe+i 



= u 



fc+1 
■U+i)' 



P 



fc+i 




If the algorithm converges, then Hindoo uj^ 1 = u fc+1 and linij-xjo WX 1 — ?7 fc+1 , 
where u fc+1 and rj k+1 are the approximate solution of the coupled problem at time step 



The algorithm entails, at each sub-iteration, the computation of the equation for the 
structure mechanics, the Navier-Stokes equations and the solution of two Laplace equations 
(12.81 . one for every displacement component. It is therefore quite computationally expen- 
sive. Alternatively, less implicit formulations may be adopted, see for instance 1131 . yet it 
has been found that for the problem at hand a strong coupling between fluid and structure 
must be maintained also at discrete level in order to have stable algorithms 1 12 1. 

3. Multiscale modelling of the cardiovascular system 

The cardiovascular system is highly integrated. In many cases, to isolate the part 
of interest from the rest of the system would require specification of point-wise boundary 
data on artificial boundary sections. These are difficult to pre-determine. To account for 
the effect of the global circulatory system when focusing on specific regions we propose 
to integrate a hierarchy of models operating at different "scales". At the highest level 
we have the full three-dimensional fluid-structure interaction problem. This will be used 
where details of local flow fields are needed. At the lowest level, we use lumped parameter 
models based on the resolution of systems of non-linear ordinary algebraic-differential 
equations for averaged mass flow and pressure. The latter models are often described by 
help of an analogy with an electrical circuit, where the voltage represents blood pressure 
and the current the flow rate. They are well fit to supply the more sophisticated models 
with the effects of the circulation in small vessels, the capillary bed, the venous system, 
as well as the action of the heart. A transition between the two extrema could be achieved 
by convenient one-dimensional models expressed by a first order non-linear hyperbolic 
system (see Fig. @}. The derivation of the one-dimensional model and a possible numerical 
implementation may be found in 16 I4l ll4l . 

An analysis of the coupling between fluid-structure models and one dimensional 
models may be found in |3|, while the direct coupling by lumped parameter models and 
Navier-Stokes model is found in 11511171 . the coupling between lumped parameter models 
and one-dimensional models is also treated in |7 1. 

Acknowledgements. The research activity here described has been supported by various 
Swiss and Italian research agencies and isntitutions, in particular the Swiss National Re- 
search Fund (FNS), and the Italian CNR, Ministry of Education (MIUR). The author also 
thanks Luca Formaggia for his contribution to the preparation of this paper. 



t k+i 



Mathematical Modelling and of the Cardiovascular System 



847 



Pi 




Left Right Pnlmunaiy 

Ventricle Ventricle System 



Figure 4: An example of multiscale simulation of blood flow, with 
the interplay between three-dimensional, one-dimensional and lumped 
parameters models. On top we show a global model of the circulatory 
system where a coronary by-pass is being simulated by a Navier-Stokes 
fluid-structure interaction model. The rest of the circulatory system is 
described by means of a lumped parameter model, based on the solution 
of a system of ODEs, is here represented by an electrical circuit analog 
in the bottom part of the figure. 
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